Departamento de Ecologia IB-USP

Um problema

Qual é a alturas média das årvores em uma floresta?

InferĂȘncia estatĂ­stica

  • Fazer afirmaçÔes sobre as caracterĂ­sticas de uma população, com base em informaçÔes dadas por amostras.

  • População Ă© qualquer conjunto de elementos que vocĂȘ estĂĄ observando

  • Amostra Ă© qualquer subconjunto desta população.

Ajustando um modelo

  • Estimar parĂąmetros da população baseados na amostra, modelo teĂłrico de distribuição de probabilidades

Modelo de distribuição para altura das årvores

  • VariĂĄvel AleatĂłria de distribuição de probabilidades Normal com mĂ©dia (ÎŒ) e desvio padrĂŁo (σ).

## [1] 22.72789
## [1] 32.3065
## [1] 5.683881

Definindo variĂĄveis aleatĂłrias

A variĂĄvel aleatĂłria X Ă© uma variĂĄvel que tem um valor Ășnico (determinado aleatoriamente) para cada resultado de um experimento (lembre-se esse termo veio da teoria de probabilidade).

Exemplos:

  1. Altura das ĂĄrvores de uma floresta (!!)

  2. NĂșmero de presas capturadas por um predador em um determinado dia;

  3. Comprimento de um peixe adulto selecionado aleatoriamente.

As variĂĄveis aleatĂłrias podem ser discretas ou contĂ­nuas.

Função de probabilidade

Função que define as probabilidades P(X) para cada valor possível da variåvel X.

Distribuição de Probabilidades.

Algumas DistribuiçÔes de Probabilidade

Bernouli: X ~ Bernouli(p)

  • Discreta, valores possĂ­veis 0 e 1. p - probabilidade de sucesso
  • Jogar moedas: cara - 0, coroa - 1.
  • Presença/ausĂȘncia de uma espĂ©cie em locais amostrados

Binomial: X ~ Bin(n, p)

  • Discreta: nĂșmero de sucessos em n tentativas . p - prob sucesso.

  • NĂșmero sementes predadas em cada experimento contendo 20 semente.

  • NĂșmero de animais infectados em relação ao nĂșmero de capturados em cada local

Exemplo Binomial

  • Predação de girinos por odonatas, p = 0.30.
    Total girinos: 6 num lago.
  • Qual probabilidades de que 0, 1, 2, 3, 5 ou todos sejam predados?

Poisson: X ~ Pois(λ)

  • Discreta: Probabilidade de uma sĂ©rie de eventos ocorrem em um perĂ­odo fixo de tempo, ĂĄrea, volume, quadrante, etc.

  • Dados de contagem

  • AbundĂąncia de espĂ©cies em um local

  • NĂșmero de indivĂ­duos capturados por tempo

Exemplo Poisson

  • NĂșmero de visitas de borboletas em uma planta em 15 min: 10
  • Probabilidade de 5 visitas dpois(5,10)

Normal: X ~ N(ÎŒ, σ)

  • SimĂ©trica, contĂ­nua

  • Metade (50%) dos dados estĂĄ acima (e abaixo) da mĂ©dia

  • 68% dados estĂĄ dentro de 1 desvio padrĂŁo da mĂ©dia

  • 95% estĂĄ dentro de 2 desvios padrĂ”es da mĂ©dia

  • Virtualmente todos os valores estĂŁo dentro de 3 desvios padrĂ”es da mĂ©dia

Exemplo Normal

  • Qual Ă© a probabilidade de que um peixe capturado aleatoriamente tenha 20,15 cm ou mais? MĂ©dia da população Ă© 17,1 cm, desvio padrĂŁo 1,21 cm?
## [1] 0.005856729

Teorema do Limite Central

"Toda soma de variĂĄveis aleatĂłrias independentes de mĂ©dia finita e variĂąncia limitada Ă© aproximadamente Normal, desde que o nĂșmero de termos da soma seja suficientemente grande."

Independentemente do tipo de distribuição da população, na medida em que o tamanho da amostra aumenta, a distribuição das médias amostrais tende a uma distribuição Normal.

Exemplo TLC

Voltando ao nosso problema

  • Dados de duas espĂ©cies
levels(arv$sp)
## [1] "sp1" "sp2"

Serå que a altura média de uma espécie nessa floresta é diferente da altura média da outra espécie?

Ou seja, serå que as espécies possuem uma mesma distribuição de probabilidades e as diferenças encontradas são devido ao acaso ou serå que cada espécie é uma variåvel aleatória com médias diferentes?

DistribuiçÔes de frequĂȘncias das alturas para cada espĂ©cie

Diferença entre as médias das alturas

mean(arv$alt[arv$sp=="sp2"]) - mean(arv$alt[arv$sp=="sp1"])
## [1] 5.109466

Posso afirmar que as médias são diferentes?

H0 - não hå diferença entre as alturas das duas espécies. A média da altura da espécie 1 é igual à média da altura da espécie 2 e a diferença observada se deve ao acaso.

H1 - hå diferença entre as alturas das duas espécies. A média da altura da espécie 1 é diferente da média da altura da espécie 2.

Erros associados a escolher cada hipĂłtese como verdadeira

AFIRMAÇÃO: nĂŁo hĂĄ diferença entre as alturas das espĂ©cies,
VERDADE: hå diferenças.

AFIRMAÇÃO: hĂĄ diferença entre as alturas das espĂ©cies,
VERDADE: não hå diferenças.

Definição formal dos erros:

ERRO DO TIPO I (α): rejeitar a hipótese nula dada que ela é verdadeira.

ERRO DO TIPO II (B): aceitar a hipĂłtese nula dado que ela Ă© falsa.

Testando as hipĂłteses

Contrastar a probabilidade de que a diferença que vocĂȘ encontrou nas mĂ©dias das alturas das duas espĂ©cies pode se considerada devido ao acaso ou nĂŁo.

Qual é a probabilidade de que as duas amostras pertencem à mesma população de medidas?

Podemos ou não rejeitar a hipótese nula (de que a diferença é ao acaso), baseado em alguma probabilidade de estarmos cometendo algum erro (α).

Graus de liberdade

"O nĂșmero de observaçÔes independentes menos o nĂșmero de parĂąmetros estimados."

Quantas observaçÔes independentes foram usadas para calcular a média e a variùncia dos dados sob a hipótese nula?

Nosso exemplo: 100 ĂĄrvores - 2 parĂąmetros = 98 graus de liberdade

EstatĂ­stica t

Não entraremos em detalhes de como se calcula a estatística t e como é a sua distribuição de probabilidades (recomendo ver na wikipedia). O importante a saber aqui é que tendo as duas amostras modeladas como uma distribuição normal, calculamos a estatística t e com um certo valor de graus de liberdade, nós podemos recorrer às tabelas da distribuição t para ver qual a probabilidade de se obter aquele valor que obtivemos dada que a hipótese nula é verdadeira. Se esta probabilidade for muito pequena, a chance de estarmos caindo no erro do tipo 1 é pequena.

t.test(arv$alt~arv$sp, var.equal=T)
## 
##  Two Sample t-test
## 
## data:  arv$alt by arv$sp
## t = -5.0125, df = 98, p-value = 2.387e-06
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -7.132311 -3.086621
## sample estimates:
## mean in group sp1 mean in group sp2 
##          20.17316          25.28262
#o argumento var.equa=T diz que as variùncias das duas espécies são as mesmas
# podemos fazer o teste com a premissa de que as variĂąncias nĂŁo sejam iguais,
# usando var.equal=F, isso vai "consumir" graus de liberdade.

Modelos lineares

E se tivéssemos na verdade 3 espécies de årvores para compararmos as alturas? Acabamos de ver que o teste t serve apenas para comparar duas amostras. O que chamamos de modelos lineares é uma grande família de modelos estatísticos que modelam uma relação linear entre a variåvel dependente (Y) e a(s) variåvel(is) independente(s). Nesta família encontra-se o teste t que acabamos de ver, a Anålise de Variùnci (ANOVA), a regressão, etc.

Chamamos de Modelos Lineare Gerais (LMs), aqueles que possuem como premissa a variåvel dependente Y como de distribuição Normal, e de Modelos Lineares Genearlizados (GLMs) os que assumem outros tipos de distribuição para a variåvel dependente (como a Binomial e Poisson). No R, a função que usamos para estes modelos são lm e glm (hå alguns mais específicos como o prórpio t.test e o aov para ANOVA). O R lida com todos os modelos lineares de maneira similar, e um detalhe importante é que todas as anålises desta família devem ser salvas em um objeto, para seu resultado ser apresentado após o comando summary.

É importante lembrar que neste roteiro estamos lidando com modelos paramĂ©tricos, ou seja, que modela a variĂĄvel de interesse como uma variĂĄvel aleatĂłria pertencente a uma certa distribuição de probabilidades. Existem outras abordagens estatĂ­sticas nĂŁo-paramĂ©tricas que fazem testes que nĂŁo tem como premissa a distribuição de probabilidades. Veja um exemplo de teste t nĂŁo paramĂ©trico neste roteiro.

Mesmo fazendo parte de uma mesma família de modelos, nas sessÔes seguintes vamos separar os modelos lineares pelos nomes tradicionais das anålises: anålise de variùncia, regressão linear simples e anålise de covariùncia. Depois, veremos alguns exemplos de GLMs aplicåveis a dados ecológicos.

ANOVA

Vamos entĂŁo incluir no nosso exemplo inicial uma terceira espĂ©cie amostrada (e vamos diminuir o nĂșmero de amostras por espĂ©cie para ficar mais fĂĄcil de fazer os cĂĄlculos passo a passo):

alt.sp3 <- c(36.90076761,33.23131727,32.82767311,24.93410577,31.87626329,
         28.76824248,26.6436144,31.41238398,27.65383929,31.70425719,
         33.60752445,28.30263354,29.32210674,23.4790054,22.20793046,
         18.28797293,23.39044736,28.75820917,29.94344703,34.47430235,
         32.09820862,26.99845592,27.64858951,31.90557306,33.63236368,
         25.35966179,35.00885172,32.89125598,18.98802229,30.8516906,
         22.04191398,35.37232889,35.67150408,34.30257037,24.43038475,
         20.6274663,31.15717916,23.27023231,36.20568545,34.33068448,
         28.40769054,31.00578874,29.50989158,31.69915991,37.25975797,
         30.1279997,38.623371,35.36472829,29.88304259,31.00145491)

sp3 <- data.frame(alt=alt.sp3, sp=rep("sp3",50))
arv2 <- rbind(arv,sp3)

amostras <- seq(1,150,by=5)

arv3 <- arv2[ amostras, ]

Neste exemplo lidaremos com um modelo simples de anova chamado one way ANOVA, porque lida apenas com uma variĂĄvel independente categĂłrica (fator).

A ANOVA testarå a hipótese nula de que as médias das alturas das 3 espécies não diferem.

Faremos agora a ANOVA "na unha", para entendermos melhor cada passo da anålise e a interpretação de seus resultados para depois usarmos a função do R que faz esta anålise diretamente.

ANOVA na unha

Vamos calcular os valores da tabela de ANOVA. Começando com os desvios quadråticos, ou seja, quanto os dados desviam da média (idéia parecida com a variùncia). O ponto importante é que essa variação é aditiva e portanto, pode se decomposta. A variação total é decomposta em variação relacionada ao tratamento (entre grupos), no nosso caso às espécies, e uma variação interna dos grupos (chamada de erro). A estatística F é a razão entre essas variaÔes após dividir cada uma delas pelos seus respectivos graus de liberdade. Complicou? Vamos fazer os cålculos e ver os gråficos para ver se entendemos melhor:

A tabela de ANOVA que vamos construir Ă© essa:

Primeiro vamos mudar a forma do nosso data frame para facilitar os comandos

arv4 <- data.frame(sp1=arv3$alt[arv3$sp=="sp1"],
                   sp2=arv3$alt[arv3$sp=="sp2"],
                   sp3=arv3$alt[arv3$sp=="sp3"])
arv4
##          sp1      sp2      sp3
## 1  21.282000 19.14508 36.90077
## 2  19.875790 18.66252 28.76824
## 3  23.562570 25.10581 33.60752
## 4  19.313744 18.10266 18.28797
## 5  26.308616 21.48084 32.09821
## 6   8.301795 28.08559 25.35966
## 7  27.192502 20.47805 22.04191
## 8  15.941866 29.22927 20.62747
## 9  20.596624 19.32911 28.40769
## 10 22.048945 18.21446 30.12800
boxplot(arv4)

Vamos calcular a média geral, as médias para cada espécie, e as diferenças entre a média geral e a média para cada espécie, para chegar ao valor da soma dos quadrados total:

media.geral <- mean(arv3$alt)
media.geral
## [1] 23.28284
medias.sps <- apply(arv4, 2, mean)
medias.sps
##      sp1      sp2      sp3 
## 20.44245 21.78334 27.62274
dif.geral <- arv4 - media.geral
ss.especies <- dif.geral^2
ss.especies
##                sp1       sp2        sp3
##  [1,]   4.00337202 17.121094 185.447876
##  [2,]  11.60801175 21.347420  30.089610
##  [3,]   0.07824755  3.323227 106.599051
##  [4,]  15.75374565 26.834280  24.948725
##  [5,]   9.15530094  3.247214  77.710675
##  [6,] 224.43179074 23.066347   4.313177
##  [7,]  15.28543337  7.866884   1.539904
##  [8,]  53.88994194 35.359990   7.051024
##  [9,]   7.21577099 15.631997  26.264064
## [10,]   1.52250306 25.688498  46.856173
ss.total <- sum(ss.especies)
ss.total
## [1] 1033.251

Para calcular os desvios quadråticos totais, nós subtraímos cada altura das årvores pela média geral e elevamos ao quadrado. Veja o gråfico:

vetor.cor <- rep(1:3, each=10) #vetor de cores

plot(x = c(1:30), y = arv3$alt, ylim=c(10,40),
     pch=(rep(c(15,16,17), each=10)),
     col=vetor.cor,
     ylab="Variåvel Resposta", xlab="ObservaçÔes",
     main="Variação total")
    for(i in 1:30)
    {
    lines(c(i,i),c(arv3$alt[i],mean(arv3$alt)),col=vetor.cor[i])
    }
    abline(h=media.geral)

Agora vamos fazer a somatĂłria dos desvios quadrĂĄticos dentro de cada grupo (ss.intra).

O grĂĄfico para entender esse cĂĄlculo:

vetor.medias<-rep(medias.sps, each=10)

plot(c(1:30), arv3$alt, ylim=c(10,40),
     pch=(rep(c(15,16,17),each=10)),
     col=vetor.cor,
     main="Variação Intra Grupos",
     ylab="Variåvel Resposta", xlab="ObservaçÔes")
    for(i in 1:30)
    {
    lines(c(i,i),c(vetor.medias[i],arv3$alt[i]),col=vetor.cor[i])
    }
    lines(c(1,10),c(medias.sps[1],medias.sps[1]),col=1)
    lines(c(11,20),c(medias.sps[2],medias.sps[2]),col=2)
    lines(c(21,30),c(medias.sps[3],medias.sps[3]),col=3)

CĂĄlculo dos valores:

medias.sps
##      sp1      sp2      sp3 
## 20.44245 21.78334 27.62274
ss.sp1=sum((arv4$sp1-medias.sps["sp1"])^2)
ss.sp1
## [1] 262.2655
ss.sp2=sum((arv4$sp2-medias.sps["sp2"])^2)
ss.sp2
## [1] 157.0018
ss.sp3=sum((arv4$sp3-medias.sps["sp3"])^2)
ss.sp3
## [1] 322.4728
ss.intra=ss.sp1+ss.sp2+ss.sp3
ss.intra
## [1] 741.7401

A soma dos quadrados entre grupos:

plot(c(1:30), vetor.medias, ylim=c(10,40), 
  pch=(rep(c(15,16,17),each=10)),
  col=vetor.cor, 
  main="Variação Entre Grupos", 
  ylab="Variåvel Resposta", xlab="ObservaçÔes")
 for(i in 1:30)
    {
    lines(c(i,i),c(vetor.medias[i],mean(vetor.medias)),col=vetor.cor[i])
    }
    abline(h=media.geral)
    points(c(1:30),arv3$alt, ylim=c(10,50), 
           pch=(rep(c(0,1,2),each=10)), col=vetor.cor, cex=0.5)

medias.sps
##      sp1      sp2      sp3 
## 20.44245 21.78334 27.62274
media.geral
## [1] 23.28284
ss.entre=10*sum((medias.sps-media.geral)^2)
ss.entre
## [1] 291.5112
#conferindo os cĂĄlculos
ss.intra+ss.entre
## [1] 1033.251
ss.total
## [1] 1033.251

CĂĄlculo do F

# Desvios médios
ms.entre=ss.entre/2
ms.intra=ss.intra/27
ms.entre
## [1] 145.7556
ms.intra
## [1] 27.47186
#F - razĂŁo das variĂąncias
F.sps=ms.entre/ms.intra
F.sps
## [1] 5.305634

Distribuição F

Vamos ver na distribuição F qual a probabilidade de termos econtrado o valor F.sps ao acaso:

curve(expr=df(x, 2,27),
      main="Distribuição F de Fisher (df=2,27)", 
      xlab="Valor F",
      ylab="Densidade ProbabilĂ­stica (df)",
      xlim=c(2,12))
abline(v=F.sps, col="red")
abline(h=0, lty=2)

xf=seq(F.sps, 12, 0.01)
ydf=df(xf, 2, 27)
polygon(c(F.sps,xf),c(0,ydf),col="red")

text(x= 7,y=0.08,paste("pf(x) =",
     round(pf(F.sps,2,27,lower.tail=F),4)), 
     cex=0.8, col="red")

CĂĄlculo do P

p.sps=pf(F.sps, 2, 27, lower.tail=FALSE)
p.sps
## [1] 0.01139248

Para mais informaçÔes sobre o F e as comparaçÔes com o teste t, veja esse roteiro.

A tabela final

Colocando os dados calculados em nossa tabela de anova

ANOVA no R

A função que constrói esta tabela de ANOVA e faz o teste estatístico é a aov. Vamos comparar:

anova.sps <- aov(alt~sp, data=arv3)
summary(anova.sps)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## sp           2  291.5  145.76   5.306 0.0114 *
## Residuals   27  741.7   27.47                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Agora nĂłs sabemos de ondem vem os nĂșmeros nesta tabela e como foi calculado o P, certo?

Testes Ă  posteriori

Onde estão as diferenças?

  • veja os grĂĄficos

  • faça teste Ă  posteriori - Tukey HDS

TukeyHSD(anova.sps)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = alt ~ sp, data = arv3)
## 
## $sp
##             diff        lwr       upr     p adj
## sp2-sp1 1.340893 -4.4708805  7.152667 0.8360276
## sp3-sp1 7.180300  1.3685259 12.992073 0.0132022
## sp3-sp2 5.839406  0.0276327 11.651180 0.0487468
#revendo grĂĄfico dos dados
stripchart( arv3$alt~arv3$sp,
           vertical = TRUE, pch = 16, 
           col = c("black", "red", "green"))

#boxplot
boxplot(arv4)

Outra maneira de fazer a ANOVA no R:

lm.anova.sp <- lm(alt~sp, data=arv3)
summary.aov(lm.anova.sp) #mesmo que o summary do aov
##             Df Sum Sq Mean Sq F value Pr(>F)  
## sp           2  291.5  145.76   5.306 0.0114 *
## Residuals   27  741.7   27.47                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# mostra o resultado de forma um pouco diferente, consegue entender?
summary(lm.anova.sp) 
## 
## Call:
## lm(formula = alt ~ sp, data = arv3)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.1407  -3.0002  -0.0742   3.2719   9.2780 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   20.442      1.657  12.334 1.32e-12 ***
## spsp2          1.341      2.344   0.572  0.57202    
## spsp3          7.180      2.344   3.063  0.00492 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.241 on 27 degrees of freedom
## Multiple R-squared:  0.2821, Adjusted R-squared:  0.229 
## F-statistic: 5.306 on 2 and 27 DF,  p-value: 0.01139

RegressĂŁo linear simples

  • Y contĂ­nuo - X continuo

Modelo da RegressĂŁo

alt.nutr <- lm(alt ~ nutri, data = arv3)
summary(alt.nutr)
## 
## Call:
## lm(formula = alt ~ nutri, data = arv3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.7176 -1.7171  0.0088  1.4124  9.2072 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 14.38711    1.05832  13.594 7.44e-14 ***
## nutri        0.58006    0.05971   9.714 1.82e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.906 on 28 degrees of freedom
## Multiple R-squared:  0.7712, Adjusted R-squared:  0.763 
## F-statistic: 94.37 on 1 and 28 DF,  p-value: 1.818e-10

Plot da RegressĂŁo

AnĂĄlise de CovariĂąncia

  • Y continuo

  • X1 contĂ­nuo - X2 categĂłrico

Modelo da ANCOVA: aditivo

alt.sp.nutr <- lm(alt ~ nutri + sp, data = arv3)
summary(alt.sp.nutr)
## 
## Call:
## lm(formula = alt ~ nutri + sp, data = arv3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.6682 -0.9554 -0.1585  1.3581  7.6134 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  13.3971     1.1034  12.142 3.23e-12 ***
## nutri         0.5256     0.0561   9.369 8.09e-10 ***
## spsp2         1.6421     1.1423   1.437  0.16251    
## spsp3         3.8326     1.1965   3.203  0.00357 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.553 on 26 degrees of freedom
## Multiple R-squared:  0.836,  Adjusted R-squared:  0.817 
## F-statistic: 44.16 on 3 and 26 DF,  p-value: 2.401e-10

GrĂĄfico ANCOVA aditivo

Modelo ANCOVA com interação

alt.sp.nutr3 <- lm(alt ~ nutri*sp, data=arv3)
summary(alt.sp.nutr3)
## 
## Call:
## lm(formula = alt ~ nutri * sp, data = arv3)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.6724 -0.9444 -0.3313  1.0416  7.0702 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 13.40167    1.42727   9.390 1.66e-09 ***
## nutri        0.52527    0.08906   5.898 4.38e-06 ***
## spsp2        2.93717    1.97080   1.490    0.149    
## spsp3        0.43636    2.75666   0.158    0.876    
## nutri:spsp2 -0.10095    0.12423  -0.813    0.424    
## nutri:spsp3  0.17187    0.14350   1.198    0.243    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.474 on 24 degrees of freedom
## Multiple R-squared:  0.8578, Adjusted R-squared:  0.8282 
## F-statistic: 28.96 on 5 and 24 DF,  p-value: 2.003e-09

Gråfico ANCOVA com interação

Tabela de ANOVA do modelo

summary.aov(alt.sp.nutr3)
##             Df Sum Sq Mean Sq F value   Pr(>F)    
## nutri        1  796.8   796.8 130.179 3.53e-11 ***
## sp           2   66.9    33.5   5.467   0.0111 *  
## nutri:sp     2   22.6    11.3   1.846   0.1796    
## Residuals   24  146.9     6.1                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Modelos lineares generalizados

  • Outras distribuiçÔes de probabilidades alĂ©m da Normal

  • VariĂąncias nĂŁo homogĂȘneas

Mais comuns:

  • Poisson – Ășteis para dados de contagem

  • Binomial – Ășteis para dados com proporçÔes, ou presença/ausĂȘncia

Etapas do GLM

  1. Uma hipótese sobre a distribuição da variåvel resposta Yi.

  2. Especificação do modelo (variåveis X).

3.Função de ligação entre a média Yi e as variåveis X. Função que lineariza a regressão

FunçÔes de ligação

Distribuição link function
normal identity
poisson log
binomial logit
Gamma reciprocal

GLM Poisson

  • Y - NĂșmero de atropelamentos de anuros em uma estrada

  • X - distĂąncia a um Parque Nacional.

Modelo GLM Poisson

mod <- glm(TOT.N ~ D.PARK, data = atrop, family = poisson)
summary(mod)
## 
## Call:
## glm(formula = TOT.N ~ D.PARK, family = poisson, data = atrop)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -8.1100  -1.6950  -0.4708   1.4206   7.3337  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.316e+00  4.322e-02   99.87   <2e-16 ***
## D.PARK      -1.059e-04  4.387e-06  -24.13   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1071.4  on 51  degrees of freedom
## Residual deviance:  390.9  on 50  degrees of freedom
## AIC: 634.29
## 
## Number of Fisher Scoring iterations: 4

Testando contra hipĂłtese nula

anova(mod, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: poisson, link: log
## 
## Response: TOT.N
## 
## Terms added sequentially (first to last)
## 
## 
##        Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                      51     1071.4              
## D.PARK  1   680.55        50      390.9 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

GrĂĄfico Modelo GLM Poisson

##   (Intercept)        D.PARK 
##  4.3164849800 -0.0001058506

GLM Binomial

  • Y - Quantidade de veados infectados dentre os capturados

  • X - quantidade de vegetação aberta

Modelo GLM binomial

nao.infectado <- veados$amostrado- veados$infectado
y <- cbind(veados$infectado, nao.infectado)
mod.veado <- glm(y ~  openland, data=veados, family="binomial")
summary(mod.veado)
## 
## Call:
## glm(formula = y ~ openland, family = "binomial", data = veados)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.8052  -1.1463   0.6707   1.3416   3.5483  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   2.4820     0.1736   14.29   <2e-16 ***
## openland     -6.0336     0.4432  -13.61   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 376.18  on 31  degrees of freedom
## Residual deviance: 121.30  on 30  degrees of freedom
## AIC: 204.65
## 
## Number of Fisher Scoring iterations: 5

Teste GLM binomial

anova(mod.veado, test = "Chisq")
## Analysis of Deviance Table
## 
## Model: binomial, link: logit
## 
## Response: y
## 
## Terms added sequentially (first to last)
## 
## 
##          Df Deviance Resid. Df Resid. Dev  Pr(>Chi)    
## NULL                        31     376.18              
## openland  1   254.88        30     121.30 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

GrĂĄfico GLM binomial

# salvando os coeficientes
coefs <- coef(mod.veado)
coefs
## (Intercept)    openland 
##    2.482035   -6.033583
#lembre que estes coefs foram linearizados com a função de ligação
# para a binomial a link function Ă© a logit
# então temos que transformar os y no inverso do logit, a função é essa:
invlogit = function(x) return(exp(x)/(1+exp(x)))

Checando as premissas dos modelos

Teste de normalidade

Teste de Shapiro-Wilks

  • H0: dados sĂŁo normais
shapiro.test(arv3$alt)
## 
##  Shapiro-Wilk normality test
## 
## data:  arv3$alt
## W = 0.96262, p-value = 0.3609
shapiro.test(arv$alt[arv$sp == "sp1"])
## 
##  Shapiro-Wilk normality test
## 
## data:  arv$alt[arv$sp == "sp1"]
## W = 0.9906, p-value = 0.9591
shapiro.test(arv$alt[arv$sp == "sp2"])
## 
##  Shapiro-Wilk normality test
## 
## data:  arv$alt[arv$sp == "sp2"]
## W = 0.97453, p-value = 0.3502

Teste de homogeneidade de variĂąncias

Teste de Bartlett:

bartlett.test(arv3$alt~arv3$sp)
## 
##  Bartlett test of homogeneity of variances
## 
## data:  arv3$alt by arv3$sp
## Bartlett's K-squared = 1.1109, df = 2, p-value = 0.5738

Teste homogeneidade de variĂąncias

Teste de Levene:

leveneTest(arv3$alt~arv3$sp)
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2   0.512  0.605
##       27

Checagem do ajuste do modelo - ResĂ­duos

RESÍDUO = Y(observado) - Y(dados ajustados)

A premissa de normalidade na verdade recai sobre os resĂ­duos do modelo.

resid(anova.sps)
##           1           6          11          16          21          26 
##   0.8395548  -0.5666556   3.1201253  -1.1287012   5.8661704 -12.1406501 
##          31          36          41          46          51          56 
##   6.7500566  -4.5005793   0.1541789   1.6065001  -2.6382600  -3.1208224 
##          61          66          71          76          81          86 
##   3.3224765  -3.6806771  -0.3024985   6.3022481  -1.3052922   7.4459310 
##          91          96         101         106         111         116 
##  -2.4542276  -3.5688778   9.2780228   1.1454976   5.9847796  -9.3347719 
##         121         126         131         136         141         146 
##   4.4754638  -2.2630830  -5.5808309  -6.9952785   0.7849457   2.5052549
#residuals(anova.sps) # o mesmo que a função de cima

plot do modelo

ResĂ­duos regressĂŁo

par(mfrow=c(2,2))
# regressĂŁo: alt ~ nutri
plot(alt.nutr)

ResĂ­duos ANCOVA aditivo

par(mfrow=c(2,2))
# ANCOVA modelo aditivo: alt ~ nutri + sps
plot(alt.sp.nutr)

ResĂ­duos ANCOVA interativo

par(mfrow=c(2,2))
# ANCOVA modelo interativo: alt ~ nutri * sps
plot(alt.sp.nutr2)

ResĂ­duos GLM poisson

par(mfrow=c(2,2))
# GLM poisson: atrop ~ dist.park
plot(mod)

ResĂ­duos LM com a poisson

par(mfrow=c(2,2))
plot(lm(TOT.N ~ D.PARK, data = atrop))

ResĂ­duos GLM binomial

par(mfrow=c(2,2))
# GLM binoimial: prop.infectado ~ openland
plot( mod.veado)

ResĂ­duos LM com a binomial

par(mfrow=c(2,2))
plot( lm(prop.veados ~ veados$openland))